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Collisional evolution of eccentric planetesimal swarms 



o 
o 

(N 
-i— > 

o 
O 
in 

(N 

w 
43 

9 L1, 
6 

03 



> 

in 

(N 

d 

On 
O 



13 



M. C. Wyatt 1 *, M. Booth 1 , M. J. Payne 1 , L. J. Churcher 1 

1 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CBS 0HA, UK 



25 October 2009 



ABSTRACT 

Models for the steady state collisional evolution of low eccentricity planetesimal belts identify 
debris disks with hot dust at 1AU, like rj Corvi and HD69830, as anomalous since collisional 
processing should have removed most of the planetesimal mass over their > 1 Gyr lifetimes. 
This paper looks at the effect of large planetesimal eccentricities (e 3> 0.3) on their col- 
lisional lifetime and the amount of mass that can remain at late times Mi atc . Assuming an 
axisymmetric planetesimal disk with common pericentre distances and eccentricities e, we 
find that Mi atc oc e _5 / 3 (l + e) 4 / 3 (l — e) -3 . For a scattered disk-like population (i.e., with 
common pericentre distances but range of eccentricities), in the absence of dynamical evo- 
lution, the mass evolution at late times would be as if only planetesimals with the largest 
eccentricity were present in the disk. Despite the increased remaining mass, higher eccentric- 
ities do not increase the amount of hot emission from the collisional cascade until e > 0.99, 
partly because most collisions occur near pericentre thus increasing the dust blow-out diam- 
eter. However, at high eccentricities (e > 0.97) the blow-out population extending outwards 
from pericentre may be detectable above the collisional cascade; higher eccentricities also 
increase the probability of witnessing a recent collision. All of the imaging and spectroscopic 
constraints for rj Corvi can be explained with a single planetesimal population with pericentre 
at 0.75AU, apocentre at 150AU, and mass 5M®; however, the origin of such a high eccentric- 
ity population remains challenging. The mid-infrared excess to HD69830 can be explained 
by the ongoing destruction of a debris belt produced in a recent collision in an eccentric plan- 
etesimal belt, but the lack of far-infrared emission would require small bound grains to be 
absent from the parent planetesimal belt, possibly due to sublimation. The model presented 
here is applicable wherever non-negligible planetesimal eccentricities are implicated and can 
be readily incorporated into N-body simulations. 

Key words: circumstellar matter - stars: planetary systems: formation. 



1 INTRODUCTION 

A natural byproduct of the planet formation process, at least in 
the core accretion model, is the formation of planetesimals (Lis- 
sauer 1993). Evidence for planetesimals following the protoplane- 
tary disk phase comes from debris disks, a phenomenon in which 
main sequence stars exhibit an infrared excess which is attributed 
to the thermal emission of dust released from planetesimals in col- 
lisions and sublimation (see review in Wyatt 2008). The Solar Sys- 
tem has its own debris disk, the majority of which is in the asteroid 
and Kuiper belts. 

Typically extrasolar debris disks are observed to lie in a ring 
at a single radius (Greaves et al. 2005; Kalas, Graham & Clampin 
2005; Schneider et al. 2009), or where they are not imaged the 
emission spectrum is dominated by a single temperature (Chen et 
al. 2006). This motivates considering these disks as planetesimal 
belts that are directly analogous to the asteroid and Kuiper belts 
(Moro-Martm et al. 2008), and the disks where dust is detected at 
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multiple radii (Wyatt et al. 2005; Absil et al. 2006; Smith et al. 
2009a; Backman et al. 2009; Chen et al. 2009) are usually inferred 
to have multiple planetesimal belts. In the absence of other dynam- 
ical processes, the evolution of these belts is expected to be domi- 
nated by collisions which grind away the mass of the largest objects 
into dust which is subsequently removed by radiation pressure (or 
P-R drag in the case of the Solar System) (Wyatt 2009). 

The collisional evolution of the planetesimal belts of the So- 
lar System has been studied extensively. Collision rates can be 
derived accurately between objects moving on given orbits (Opik 
1951; Wetherill 1967; Greenberg 1982; Bottke et al. 1994; Vedder 
1996; DeirOro & Paolicchi 1998), and the steady state size dis- 
tribution of the belts resulting from their collisional evolution is 
both well understood analytically (Dohnanyi 1969; Tanaka et al. 
1996; O'Brien & Greenberg 2003; Kobayashi & Tanaka 2009) and 
one-dimensional numerical models of this evolution that include a 
realistic prescription for the outcome of collisions provide a good 
fit to the observed size distributions in the asteroid belt (Davis et 
al. 1989; Durda, Greenberg & Jedicke 1998; Bottke et al. 2005; 
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O'Brien & Greenberg 2005) and Kuiper belt (Stern & Colwell 
1997; Davis & Farinella 1997; Kenyon & Bromley 2004). 

The approach to considering the collisional evolution of extra- 
solar debris disks is slightly different in that the orbital element and 
size distributions of the parent bodies are poorly constrained, rather 
it is important to generalize the effect of this evolution on debris 
disk observability with respect to parameters such as initial plan- 
etesimal belt mass, radius and mean eccentricity. Analytical mod- 
els that achieve this were developed by Dominik & Decin (2003) 
who considered the evolution of a mono-disperse planetesimal belt 
(i.e., with planetesimals all of the same size) that feeds a population 
of smaller planetesimals and dust that is observed. This model was 
later refined by Wyatt et al. (2007a) to consider the parent planetes- 
imals and the smaller objects to form a continuous size distribu- 
tion defined by a single power law as expected for the steady state 
case where planetesimal strength is independent of size (Dohnanyi 
1969; Tanaka et al. 1996). A size dependent planetesimal strength 
was later included in such models by Lohne et al. (2008) result- 
ing in a more realistic 3 phase size distribution (e.g., O'Brien & 
Greenberg 2003). Both the Wyatt et al. (2007b) and Lohne et al. 
(2008) models were applied to the statistics of detections of debris 
disks around A stars and Sun-like stars to show that these could be 
explained if the majority of such debris disks evolve purely due to 
steady state collisional evolution. 

One important result that came out of the Wyatt et al. (2007a) 
study was the concept of a maximum planetesimal belt mass, and 
hence a maximum dust luminosity, that can remain for a given ra- 
dius belt at a given time, regardless of its initial mass. Although this 
is no longer strictly true when a size dependent strength is used, 
Lohne et al. (2008) showed that initial mass has a relatively mod- 
est effect on the mass remaining at late times, and concurred that 
for realistic planetesimal belt parameters there is indeed a maxi- 
mum planetesimal belt mass and dust luminosity for a given age 
(see also Heng & Tremaine 2009). This concept was used by Wy- 
att et al. (2007a) to show that l-2Gyr systems like r\ Corvi and 
HD69830 that have large quantities of hot dust at ~ 1 AU (Wyatt et 
al. 2005; Beichman et al. 2005; Smith, Wyatt & Dent 2008), cannot 
be replenishing that dust from planetesimal belts that are coinci- 
dent with the dust (i.e., analogous asteroid belts). They concluded 
that the parent bodies of the observed dust must have originated 
at larger radii (3> several AU) where collisional processing times 
would have been longer. The paper also concluded that the hot dust 
is transient and proposed that this might have been scattered in from 
an outer belt in an epoch akin to the Late Heavy Bombardment in 
the Solar System (see review in Hartmann et al. 2000). There are 
now several examples of systems exhibiting hot dust that appears to 
be transient by the criterion described by Wyatt et al. (2007a) (e.g., 
di Folco et al. 2007; Akeson et al. 2009; Moor et al. 2009). 

The motivation of this paper is to consider whether it is possi- 
ble to circumvent the conclusion that the hot dust in systems like r\ 
Corvi and HD69830 must be transient by postulating a population 
of parent planetesimals on highly eccentric orbits (e 3> 0.3). In 
such a model the hot dust would originate from material close to 
pericentre, and the parent population could be long-lived because 
the planetesimals spend most of their time at apocentre. This would 
challenge our traditional view of debris disks as belts of planetesi- 
mals with modest eccentricity (e < 0.3), which is also implicit in 
the models of Wyatt et al. (2007a) where collision velocities are as- 
sumed to be proportional to Keplerian velocity times a mean eccen- 
tricity for the belt, and in the models of Lohne et al. (2008) where 
eccentricities up to 0.35 were considered. However, it is clear from 
the Solar System that there are also populations of planetesimals 



on highly eccentric orbits (e ^> 0.3): the comets scattered in from 
the Kuiper belt (Duncan 2008); the scattered disk of the Kuiper 
belt (which may be primordial in origin and extends all the way 
to the Oort cloud) (Gomes et al. 2008); as well as the Near Earth 
Asteroids (Bottke et al. 2002). While the contribution of these pop- 
ulations to the dust content of the zodiacal cloud may be small, the 
cometary contribution could be as much as 90% (Ipatov et al. 2008; 
Nesvorny et al. 2009), and may have been significantly higher in the 
past, e.g., during the epoch known as the Late Heavy Bombardment 
(Gomes et al. 2005; Booth et al. 2009). Furthermore the opposite 
may be true for planetary systems with different architectures and 
formation scenarios, in which eccentric planetesimals may domi- 
nate. Indeed, planet formation models often predict a highly eccen- 
tric remnant planetesimal population (Edgar & Artymowicz 2004; 
Mandell, Raymond & Sigurdsson 2007; Payne et al. 2009). 

Thus, here we develop the model of Wyatt et al. (2007a) to 
include interactions between planetesimals of arbitrary eccentric- 
ities and semimajor axes (and inclinations). Although this model 
does not (yet) include the more realistic assumption of a size de- 
pendent planetesimal strength, it benefits by providing simple ana- 
lytical formulae for collision lifetimes from which the observability 
of a planetesimal belt as a function of its eccentricity can be read- 
ily assessed. The inclusion of a size dependent strength would be 
expected to give results within an order of magnitude of those pre- 
sented here (see, e.g., Fig. 1 1 of Lohne et al. 2008), a level of un- 
certainty that is commensurate with the uncertainty in estimates for 
planetesimal strength at each size for the Solar System and for dif- 
ferent assumptions about planetesimal composition (see e.g., Fig. 1 
of Durda et al. 1998, and Fig. 1 1 of Leinhardt & Stewart 2009). 

In ij2]we consider the collisional evolution of an axisymmet- 
ric disk of planetesimals all of which have the same pericentre and 
apocentre distances, and show that the concept of a maximum re- 
maining mass for a given age also applies in this case, but that the 
remaining mass is higher for higher eccentricities (if pericentre dis- 
tance is kept constant). The approach to calculating collision rates 
is similar to that of Bottke et al. (1994) in that we assume ran- 
dom mean longitudes, arguments of pericentre and longitudes of 
ascending node, but differs in using a particle-in-a-box approach 
to calculate the collision rate at a particular point on the orbit then 
integrating around the orbit (as opposed to calculating this from 
the fraction of the orbits that the planetesimals spend close enough 
that they overlap in physical space). Our collision rate at each point 
also includes an integration over the size distribution of impactors 
that can cause a catastrophic collision, whereas this integration is 
performed after calculating the intrinsic collision probability by 
Bottke et al. (1994), requiring that method to keep track of the 
velocity probability distribution. Our integration is performed us- 
ing a Monte-Carlo technique, but in the case of mutual collisions 
amongst a population with common eccentricities and semimajor 
axes, and assuming that collision velocities are dominated by ra- 
dial motion (due to eccentricities) rather than vertical motion (due 
to inclinations), the collision rate can also be derived analytically. 

In Sj3] we consider the more general situation in which plan- 
etesimals can interact with planetesimals with different pericentre 
and apocentre distances, and show that our collision rates agree 
with those of the most accurate studies available in the literature. 
To consider the evolution of a realistic planetesimal belt where a 
range of eccentricities and semimajor axes is present we adopt an 
approach similar to Rrivov et al. (2005, 2006) in that we consider 
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the evolution of the phase space distribution^ Here we outline a 
scheme for evolving the phase space distribution numerically, and 
apply this to a scattered-disk like population with common pericen- 
tre and a range of apocentre distances. The emission properties of 
eccentric rings are considered in Sj4]to assess if the emission spec- 
trum of real systems can be consistently explained by steady state 
processing given the stellar ages. The conclusions are given in [J5] 



2 COLLISIONAL LIFETIME OF SINGLE 
PERICENTRE-APOCENTRE POPULATION 

Our approach to calculating collision rates is based on the particle- 
in-a-box approach, wherein the planetesimals are assumed to be 
spread uniformly around an annulus and to have a mean collision 
velocity (e.g., equation 28 of Wetherill 1967). The resulting col- 
lision rate is accurate to a factor of 2 when the asteroid belt is 
considered as a single annulus (Wetherill 1967). Such techniques 
have also been well developed for studies of the accumulation of 
planetesimals into planets (e.g., Greenberg et al. 1978; Wetherill & 
Stewart 1989) and it is possible to derive precise collision rates for 
certain assumptions about the distributions of planetesimal eccen- 
tricities and inclinations, such as that these follow a Rayleigh dis- 
tribution (Greenzweig & Lissauer 1990; Lissauer & Stewart 1993). 
However the derived collision rates are only valid where eccentric- 
ities and inclinations are <C 1. Here we consider collision rates 
between pairs of orbits that can (but do not necessarily) have sig- 
nificant eccentricity by splitting up the orbits into annuli, since the 
density and velocity distributions in each annulus are well defined 
and the particle-in-a-box approach can be used to work out accu- 
rate collision rates which can then be integrated around the orbit; 
in this respect our approach to calculating collision rates is similar 
to that of Spaute et al. (1991) which assumed e,/ < 1. 

2.1 Local collision rates 

Consider a planetesimal of diameter D that is moving through a 
disk of planetesimals with a range of sizes, where the size distri- 
bution is defined such that <j(Am)GiAm is the fraction of the total 
cross-sectional area in the distribution that is in the size range Am 
to Di m + dDi m . If the local volume density of cross-sectional area 
of planetesimals (of all sizes) is <r v in AU 2 /AU 3 , and the relative 
velocity of collisions is v Ic i in m/s, then a particle-in-a-box calcu- 
lation gives the local rate of impacts from planetesimals in the size 
range Am to Am + dD im as R co i(D, Am)dAm where 

Rcai(D, Am) = 2.11 x HT 4 , f(D, D im )a v v rcl (1) 

in yr" 1 , where the constant here (and in later equations) arises from 
the choice of units for the various parameters, and 

f(D, Am) = a(Am)(l + D/ Am) 2 . (2) 

Note that due to the high relative velocities in an eccentric disk 
we have ignored gravitational focussing in this collision rate (see 

The majority of the collisions suffered by the planetesimal 
cause negligible mass loss. Although the cumulative effect of such 
cratering collisions can be important (e.g., Kobayashi & Tanaka 
2009), here we consider the rate at which the planetesimal suffers 

1 Note that the Krivov et al. 2005 model is not accurate for high eccentric- 
ities because of the way mean impact velocity was calculated (see >13.21 . 
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collisions that have sufficient energy to cause catastrophic disrup- 
tion. A catastrophic collision is defined as one in which the largest 
fragment remaining following the collision (due to both shattering 
and subsequent gravitational reaccumulation) has half the mass of 
the original planetesimal, and a planetesimal's dispersal threshold 
(in J kg -1 ) is defined such that catastrophic collisions are those 
in which the specific incident kinetic energy exceeds Qd- Disper- 
sal thresholds have been ascertained for planetesimals of varying 
size and composition using a combination of laboratory experi- 
ments and numerical modelling (e.g., Fujiwara et al. 1989; Benz & 
Asphaug 1999). This definition means that for a given collision ve- 
locity there is a minimum size of object that the planetesimal must 
be impacted by to be catastrophically destroyed. Denoting this as 
Ac gives the minimum relative size of impactor to target for catas- 
trophic disruption as X c — D tc /D where 

X c = (2Q* D /vl 1 ) 1/3 . (3) 

Working out the rate of catastrophic collisions, R cc , then requires 
integrating equation l[T) from Ac up to the largest object in the size 
distribution 

R CC (D) = 2.11 X 10- 4 / cc (AVv«rel, (4) 

where f cc (D) = ^ZiD) D ™) dD ™- 
2.2 Evolution of the size distribution 

To simplify our model for the evolution of the size distribution we 
make the assumption that the distribution follows a single power 
law 

n(D) oc D 2 ' 3qd (5) 

between sizes Anin (in /im) and Anax (in km), where the planetes- 
imals are assumed to be spherical to get <r( D) oc D 4_39d . For the 
situation where planetesimal strength (i.e., the dispersal threshold) 
is independent of size and the size distribution has no maximum or 
minimum size it is well known that the planetesimal belt's steady 
state solution has a power law size distribution with qa = 11/6 
(Dohnanyi 1969; Tanaka et al. 1996), which we use throughout this 
paper. This slope arises because it means that the mass loss rate 
in each size bin is both independent of size and balanced by mass 
gain from the fragmentation of larger planetesimals. A more real- 
istic size distribution is truncated both at small sizes (e.g., due to 
radiation pressure) and at large sizes (e.g., set by the largest object 
which formed in the belt). The truncation of the size distribution at 
the small size end causes a ripple in the steady state size distribu- 
tion (Campo-Bagatin et al. 1994). However, the truncation at large 
sizes has an important longer term effect, as it means that collisions 
eventually deplete the number of large objects, since these are no 
longer being replenished by the destruction of yet bigger objects. 

In the simple model proposed by Wyatt et al. (2007)a, the size 
distribution is considered to be in quasi-steady state thus maintain- 
ing the power law slope of equation ((5), but for mass to be lost as 
the largest planetesimals are depleted on their catastrophic collision 
timescale. For the size distribution of equation ((5), the factor in the 
collision rate for the largest planetesimals is 

/cc(D maX ) = (lO^Anin/Anaxf^G^Ac), (6) 

where the assumption that gd = 11/6 results in G(q,X c ) = 
0.2X~ 2 S + 0.67X~ 1B + X~ - 5 - 1.87 (see Wyatt et al. 2007a), 
which we further simplify to G(ll/6, X c ) » 0.2X C " 2 5 which is 
accurate to 71% for X c < 0.87 and to 30% for X c < 0.1. This 
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means that the local catastrophic collision rate of the largest plan- 
etesimals is 

flcc(Anax) M 7.49 X 10" 10 AnFn Anax' Od " ^ "oVU^ . (7) 

Since most of the objects which are causing this catastrophic 
destruction are both large and contain most of the mass (but little of 
the cross-sectional area), it is more appropriate to rewrite equation 
using the following relation 

a tot = 12650M tot p^DjfD^g (8) 

in AU 2 , where Mtot is total mass in the distribution in M®, and p 
is planetesimal density in kg m~ , to find that 

i?cc(Anax) = KM tot a v v 8 J; 3 , (9) 
where 

K = 9.5 x 10- 6 p _1 ^LQd~ 5/6 (10) 

and av = ffv/ctot is the normalised volume density of cross- 
sectional area in AU -3 . Later in the paper we adopt K = 2.1 x 
10 -14 as a fiducial value corresponding to p = 2700 kg m -3 , 
Auax = 2000 km and Qd = 200 J kg" 1 , or some other equiva- 
lent combination of these parameters. 

We acknowledge that the above prescription gives a highly 
simplistic view of the evolution of the size distribution. For ex- 
ample, for an infinite collisional cascade in which oc D a the 
size distribution is still expected to follow a power law, but with a 
slope that depends on a (O'Brien & Greenberg 2003; Kobayashi 
& Tanaka 2009). Since experiments have shown that Qjjj has a dif- 
ferent slope at small sizes (in the strength regime) to that at large 
sizes (in the gravity regime), and the primordial planetesimal distri- 
bution may differ from that expected when the distribution reaches 
steady state, a more realistic prescription for the size distribution 
has 3 power laws in different size regimes (Lohne et al. 2008), with 
numerical simulations showing that the transition from strength to 
gravity scaling also causes a further wiggle in the size distribu- 
tion (Durda et al. 1998). However, despite the seemingly simplistic 
view of the evolution presented here, it was found that its predic- 
tions for the evolution of mass and area in a planetesimal belt are 
accurate to within an order of magnitude of more detailed models 
(see, e.g., Fig. 1 1 of Lohne et al. 2008). Given that there are large 
uncertainties in the various parameters that make up the constant K 
(equation I lOt — e.g., estimates for planetesimal strength at a given 
size vary by several orders of magnitude both from constraints from 
the Solar System's debris belts and for different assumptions about 
planetesimal composition (see e.g., Fig. 1 of Durda et al. 1998, 
and Fig. 1 1 of Leinhardt & Stewart 2009) — we consider that this 
model is sufficient to assess the impact of planetesimal eccentricity 
on the observability of a debris disk, but note that it is possible to 
construct a more realistic (and more complex) prescription for the 
size distribution if assumptions about how planetesimal strength 
varies as a function of size are made. 

2.3 Low eccentricity approximation 

Collision rates calculated using particle-in-a-box methods such as 
that described above typically use what we will call the low ec- 
centricity approximation, since the assumptions break down when 
e ^> 0.3. When applied to a planetesimal belt such as the asteroid 
belt (e.g., Wetherill 1967), it is assumed that the planetesimals are 
spread uniformly throughout a torus of radius r, width dr and ver- 
tical extent 2I max r, so that normalised cross-sectional area density 
can be approximated at all locations within the torus by 



<T V = (47rr 2 dr/ ma x) \ (11) 

Second it is assumed that, although collisions in a realistic plan- 
etesimal belt occur at a range of relative velocities, collision rates 
can be calculated using a mean relative velocity; this can either be 
calculated by considering the observed distribution of orbital ele- 
ments (i.e., for the asteroid belt), or where this is not known by 
assuming that it is proportional to the Keplerian velocity at r, 

v k {r) = 2.98 x 10 4 M°- 5 r" ' 5 , (12) 

in m s , where M* is in Mq and r in AU, through the mean 
eccentricities and inclinations of planetesimals' orbits so that 

«rel/«k = /(e,I). (13) 

It can be shown that for Rayleigh distributions of eccentricities and 
inclinations, the mean relative velocities between planetesimals is 
given by f(e,I) = e^/1.25 + (I/e) 2 (equation 17 of Lissauer 
& Stewart 1993). These two assumptions can be fed directly into 
equation ((9) to get an expression for the collision rate which in 
this approximation is the same for the largest planetesimals at all 
locations within the torus: 

■Rcc(A™x) = 1.9 x 10 6 p" 1 AnLQD~ 5/6 M tot X 

r-^idr/r-r'M^I-Le^, (14) 
where it was further assumed that I = e. 

2.4 Higher eccentricities 

For higher eccentricities (e ^> 0.3) the assumptions of ^23] break 
down, because it is expected that both the cross-sectional area den- 
sity and the relative velocity of collisions are significantly different 
at different locations within the torus and so vary around the plan- 
etesimals' orbits. To calculate the collision rate between planetesi- 
mals on eccentric orbits we use an approach that differs from that 
pioneered by Opik (195 1), but show in ^3.2l that the two approaches 
get identical results. The method is based on the particle-in-a-box 
approach in that it assumes that equations {4]l and ((9} provide good 
estimates for the planetesimal's local collision rates as long as the 
volume density and relative velocity of collisions at that location 
are well known. For the orbital element distributions we consider, 
we expect these to be well characterised and to depend only on ra- 
dius (and latitude), and so the collision rate can be calculated by 
integrating the local rate along the planetesimal's orbit. 

To both work out the local collision rate and to perform this 
integration a Monte Carlo approach is used, wherein N planetesi- 
mals are chosen with orbital elements chosen randomly from given 
ranges. This section considers a disk of planetesimals all with peri- 
centre distances in the range q — dq/2 to q + dq/2 and apocentre 
distances in the range Q~dQ /2 to Q+dQ /2, where dq and dQ are 
small enough to have no effect on the result, but are kept finite for 
practical reasons and to allow easy implementation in the model of 
Sj3] The disk is assumed to be axisymmetric so that pericentre ori- 
entation vj is chosen randomly from the range to 27r, as is mean 
anomaly M. The orbital elements for these N planetesimals are 
then converted into 2-dimensional positions and velocities. 

2.4.1 2D approximation 

The 3-dimensional structure of the disk is accounted for in this 
section by assuming that vertical motion has little effect on rela- 
tive velocities which can thus be calculated in 2-dimensions, but 
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Figure 1. The factor <S2d( e ) from equation (T7J determined numerically 
using different numbers of planetesimals N for e =0.1, 0.6, 0.9, 0.96 and 
0.99, compared with that determined using N = 3 X 10 5 . For each value 
of N, the calculation was repeated 10 times to determine the lcr uncertainty 
in the value derived in this way, and this is shown by the error bars. For all 
eccentricities below 0.99 the calculation has converged for N > 3 X 10 . 
The number of radial bins was set at 40 for this calculation, which similar 
plots show is sufficient for this eccentricity range. 



that the vertical structure does affect the local cross-sectional area 
density which is accounted for by assuming a vertical extent of 
2r/ max . Given the axisymmetric structure of the disk, the collision 
rate should be the same for planetesimals at the same radius. Thus 
the disk was split into a number of annuli, iV ann , r , logarithmically 
spaced in radius between q — dq/2 and Q + dQ/2. 

The planetesimal positions were then used to derive cr(r) 
which is both the fraction of the total cross-sectional area in the 
annulus at r ± dr/2, and the fraction of time spent by planetesi- 
mals in the different annuli as they go around their orbit. Further 
dividing by the volume of the annulus, dV2d = 47rr 3 (dr/r)7 max 
gives the normalised volume density a v (r). For each planetesimal, 
its nearest neighbour is found and the difference in their velocities 
used to estimate the relative velocity of collisions. Taking relative 
velocities to the 8/3 power and averaging for each annulus results 
in (Vre\)(r). The mean collision rate for individual planetesimals 



of size D m ax in this population is then calculated using 



( )=™ ot ^a(r) f 7 v (r)K 8 / 1 3 )(r), 



(15) 



where 2~2 r ' and similar notation elsewhere, means to sum over ra- 
dius, and the extra <r(r) term (compared with equation[9} accounts 
for the fraction of time the planetesimal spends in different annuli. 

One further simplification is possible to this collision rate by 
expressing this in terms of r = r/a using the fact that a(r) = 
o"(f), and that relative velocity is proportional to the Keplerian ve- 
locity at r = a times some function of e and f (see equation[T9t so 
that 



Rcc{D : 

& d (e 



) = 



A"M tot ^ /3 (a)a' 3 (4 7 r/ max )- 1 S(e), (16) 
Y, f- 2 df- 1 a(f) 2 ([ Urcl /« k (a)] 8/3 }(f), (17) 



where S(e) = S , 2d(e) in equation l l 1 6t in the 2D approximation. 

Successful implementation of this routine requires that the 
number of planetesimals (TV) is sufficiently large for the relative 
velocity of encounters to be well approximated by the difference 
in the velocities with nearest neighbours. The number of annuli 



(-^Vann.r) must also be sufficiently large to resolve the radial vari- 
ations in collision rate between pericentre and apocentre. To as- 
sess this, Fig. [TJplots S(e) as a function of N for iVann.r = 40, 
normalised to the value expected when S(e) is calculated with 
N = 3 x 10 5 and JVann.r = 40. It is evident that the solution 
converges for TV > 10 5 for e ^ 0.99. A similar plot showing the 
effect of changing the number of radial bins shows 7V a nn,r = 40 is 
sufficient for this eccentricity range, with values as low as 10 also 
giving results within 5%. 



2.4.2 Analytical collision rates 

The 2D collision rates can also be derived analytically for a popula- 
tion with a common semimajor axis and eccentricity, as the various 
factors in equation J17| > are the consequence of 2 body Keplerian 
motion. The distribution of cross-sectional area is determined by 
the rate of radial motion a(r) jdf — 2a/(t pel dr/dt), where t per is 
the orbital period, giving: 



a/(du) — 



2 (<- 



■2u-l]~ 



(18) 



where u = a/r. The relative velocity of collisions at each radius 
has a bimodal distribution with half of collisions occurring at zero 
velocity (for planetesimals moving in the same direction), and the 
remainder at a velocity given by \/2 — 2 cos <j> times the orbital 
velocity at that radius, where (f> is the angle between the velocity 
vectors for planetesimals moving in different directions (e.g., those 
returning to pericentre colliding with those that have just passed 
pericentre), giving: 



«rcl 



/v k (a) 



l) + 2w-l]° 



(19) 



noting that {[v re i/v-k(a)] 8 )(f) is 0.5 times the expression above 
to the 8/3 power. Putting these expressions into equation (T7J and 
integrating over the range u = (1 — e) _1 to(l + e) _1 gives the 
relevant factor in the equation for the collision rate as: 



&d = 0.54e 5/3 (l-e 



-4/3 



(20) 



2.4.3 Where do collisions occur, and at what velocity? 

The first question we can answer with this model is where most 
of the collisions occur. This can be worked out from the distri- 
butions plotted in Figure l(2j. Although there is also a small den- 
sity enhancement at pericentre, planetesimals spend the majority 
of their time at apocentre (Figure [2^; see, e.g., Fig. 4b of Sykes 
1990). Nevertheless, the apocentric contribution to the overall col- 
lision rate is diminished due to both the r~ 2 term in equation ([17} 
and the higher relative velocities at lower radii (Figure [2p), where 
it should be noted that the factor in equation l |17t is relative to 
the Keplerian velocity at r — a and so is Figure |2p multiplied 
by \J a/r then to the 8/3 power. The net result is that the ma- 
jority of the collisions occur close to pericentre for all except the 
lowest eccentricities (Figure^). The Monte-Carlo simulation pro- 
vides results in excellent agreement with the analytical calculation 
( i]2.4.2t from which it can be found that 90% of collisions occur at 
r/a < (1 - e 2 )/( l - 72e). 

As noted in ^2.4.21 the distribution of collision velocities in 
this population is not uniform. Except for planetesimals that are 
close to pericentre or apocentre, relative velocities in each annulus 
have a bimodal distribution, as for planetesimals that have recently 
passed pericentre there is a low relative velocity population that 
have also recently passed pericentre and a high relative velocity 
population that are returning from apocentre to pericentre (as noted 
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Figure 2. Radial distribution of collision rate for mutual collisions amongst 
a population of planetesimals with common pericentre and apocentre dis- 
tances. The plots indicate the contributions of different factors in the colli- 
sion rate equation (T7J as a function of radius for a population with eccen- 
tricities of 0. 1 , 0.6, 0.9, 0.96, 0.99. The results of both the Monte-Carlo sim- 
ulation and the analytical calculation are plotted showing excellent agree- 
ment, (a) The distribution of cross-sectional area density, a(f)/df; the in- 
tegral under the curve for each population is 1. (b) The average collision 
velocities of planetesimals at different radii relative to the local Keplerian 
velocity [([v Icl /v k (r)] 8 / 3 )(f)] 3 / 8 . (c) The fraction of S 2 d{e) that comes 
from radii below r'. 
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Figure 3. Dependence of collision lifetime on eccentricity. For a disk in 
which mass and semimajor axis are fixed collision lifetime <x 1/S(e). The 
function S(e) is shown in (a), both calculated in the 2D approximation 
(S2d(e), solid line), and with the full 3D calculation (53d (e), dotted line), 
as well as the analytical calculation of the 2D approximation (dashed line 
which lies under the solid line). For a disk in which mass and pericentre 
distance are fixed collision lifetimes vary <x [(1 — e) 13 / 3 S(e)] — 1 . This 
function is shown in (b) with the same origin for the different linestyles as 
(a). 



in SJ2A2}. Thus the relative velocity of the high relative velocity 
population is approximately 2 3 ^ 8 times that shown in Figure |2j5, 
and remains within the range 1-2.5 times the local Keplerian ve- 
locity between apocentre and pericentre for e = 0.6 — 0.99. This 
means that, for a population with pericentre at 1AU, collisions oc- 
cur at velocities of 10 — 100 km s - , and that relatively small plan- 
etesimals can cause destructive collisions (equation[3j. 



2.4.4 Collision rate and remaining mass vs eccentricity 

Figure [3^ shows how S2d(e) varies with eccentricity, where the 
results of the Monte-Carlo simulation closely follow the predic- 
tions of the analytical calculation in equation l |20t . The factor S(e) 
can be readily used to assess how changing the eccentricity of a 
planetesimal population affects its collisional lifetime using equa- 
tion H6\ . It can also be used to consider how eccentricity affects 
the amount of disk mass that can remain at late times. For the as- 
sumptions about the evolution of the size distribution discussed in 
i)2.2| the evolution of planetesimal belt mass from its initial value 
of Mtoto can be calculated from the collisional rate using the equa- 
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tion dMtot/dt = — M to t-Rcc(-Dmax), which can be solved to give 
the mass at time t in years to be 

A/tot/M to tO = [1 + (i?cc(Anax)M tot0 /M tot )f]" 1 . (21) 

This means that at late times the remaining mass converges to a 
value of (Wyatt et al. 2007a): 

Mate = t" 1 [Mtot/i?cc(O m ax)], (22) 

= {^I^/K)[a 3 /vl ri {a)][S{e)t]-\ (23) 

where late means t S> 1/-Rcc(-D max ). 

Here we illustrate this in two ways. First we consider a disk 
in which eccentricity is varied, but the semimajor axis and disk 
mass are kept constant. The collisional lifetime of such a disk 
is oc l/5(e), meaning that increasing eccentricity results in a 
shorter collisional lifetime. The mass remaining at late times is also 
oc l/S(e), which from Fig. [3] decreases rapidly with eccentricity, 
because of the increased collision rate of material at a pericentre 
which tends to smaller radii. 

However, if pericentre location and disk mass are kept fixed 
as eccentricity is increased, collisional lifetime instead varies oc 
[(1 - e) 13/3 5(e)] _1 . Rewriting equation[23l gives 

M latc = [1.47xl0- 11 / max /^]M- 4/3 g 13/3 t- 1 [(l-e) 13/3 5(e)]- 1 : 

so that mass remaining is also oc [(1 — e) 13//3 5(e)] _1 , and this 
function is plotted in Fig. [3p. Thus both collision lifetime and 
mass remaining increases as eccentricity is increased oc e~ 5,/3 (l + 
e) 4 / 3 (l — e) -3 . Primarily this is because planetesimals spend a 
larger fraction of their orbit at large radii resulting in both a longer 
time between collisions, and a lower density of colliders in the re- 
gion where collisions occur near pericentre; as illustrated in Fig.[2] 
increasing eccentricities also results in an enhanced collision ve- 
locity at pericentre which offsets this to some extent, but there is no 
significant change in collision velocities once eccentricities have 
increased beyond 0.4. 

2.4.5 Comparison with low eccentricity approximation 

The collision rate in the low eccentricity approximation ( i]2,3t can 
also be expressed in the form of equation [T6j where a — r and 

5(e) = 5iowe(e): 

5 loW o = 3.0(dr/r) _1 e 8/3 . (25) 

This could be plotted on Fig.[3] e.g., assuming that the width of the 
ring is dr/r = 2e, which would show agreement within an order 
of magnitude for e = 0.1 — 0.98. However, this comparison is less 
instructive than noting that the masses (and fractional luminosities) 
remaining at late times in Wyatt et al. (2007a) were derived us- 
ing e = 0.05 and dr/r = 0.5 for which 5i owo — 2 x 10~ 3 and 
[(1 — e) 13 ' 3 5i owe (e)] _1 = 610. Thus to increase the amount of 
mass that can remain at late times above the values in that paper, 
assuming the radius inferred for the population corresponds to the 
location of the pericentre of an eccentric ring, would require ec- 
centricities higher than 0.78, with a factor of > 100 increase for 
e > 0.96. This illustrates the point that collisional lifetimes (and 
remaining mass) can be increased both by increasing eccentricity 
and by spreading material over a broader range of radii. 

2.4.6 3D calculation 

For the 3-dimensional calculation, it is further assumed that the 
longitudes of ascending node SI are random and that inclinations 



are randomly chosen from the range to J max , and these were 
used to calculate 3-dimensional positions and velocities for the 
planetesimals. Since collision rate is then also a function of lat- 
itude, the disk was further split in latitude into iV annj< /, bins at 
4> ± d(j>/2. The planetesimal positions were used to derive a(r, cf>), 
which is the fraction of the total cross-sectional area in the annu- 
lus at r ± dr/2 and (j> ± d(f>/2. Dividing this by the volume of the 
annulus, dV^d = 2nr 3 (dr/r) cos (f>d(f>, gives the normalised vol- 
ume density av(r, (f>). The mean velocities are also a function of <f>, 
(vfli) (r, (f))). This results in a collision rate given by equation J 1 6b 
in which 5(e) = 5 3d (e, 7 max ), where 

5 3 d(e,/ max ) = ^ r~ 2 df~ 1 y^(2J ma x/ cos (/>d(/)) x 

r <p 

a(f,0) 2 ([« rcl /« k (a)] 8/3 )(f,0). (26) 

We can now answer how important it is to account for the 
3-dimensional structure of the disk. For J max = 0.05 and using 
■Nann.^ = 9 it was found that 53d(e, I m ax)/52d(e) has a constant 
value of 1.43 ± 0.04 for eccentricities in the range 0.1 — 0.98 (see 
Figure[3}. We attribute this offset to the fact that an even distribution 
of inclinations between and J ma x results in a higher density in the 
^mid-plane (to achieve a uniform density we could have used the in- 
clination distribution function given by equation 2.20 of Krivov et 
al. 2005). Thus the 2-dimensional results give a good approxima- 
tion of the collision rates in a 3-dimensional disk, as long as the 
parameter J ma x used in this calculation is interpreted as there be- 
ing a uniform distribution of inclinations between and 1.437 m ax- 
However, for low eccentricities e <C J max , where collision veloci- 
ties are not due only to radial and azimuthal motion, but also latitu- 
dinal motion, a full 3-dimensional calculation would be necessary. 



3 EVOLUTION SCHEME FOR ECCENTRIC 
PLANETESIMAL SWARM 

The single pericentre and apocentre distance population of Sj2]was 
necessarily an idealised case given that the eccentric populations 
in the Solar System have a range of eccentricities. This section de- 
scribes how the modelling method of Sj2]can be generalised for such 
a situation, both to calculate the collision rates between populations 
with different eccentricities and semimajor axes ( ^3. It , and to use 
these rates to model the collisional evolution of a planetesimal belt 
with a distribution of orbital elements ( ^3 . 3t . The approach to mod- 
elling the collisional evolution is similar to the kinetic model em- 
ployed by Krivov et al. (2005) in that we consider the evolution 
of the phase space distribution, which here is just two dimensional 
and defined by M to t (q, Q). 



3.1 Collisions between two single pericentre-apocentre 
populations 

The parameter space is divided into cells of size q ± dq/2 and 
Q ± dQ/2, hereafter simply referred to as cell (q,Q), where 
Mtot {q, Q) is the total mass in that cell. Cells are logarithmically 
spaced. For now it is assumed that the distribution of inclinations is 
the same for each cell, and other angles are randomised as before, 
and so where needed the spatial distribution from material in each 
cell (a(r)(q, Q) and a v (r)(q, Q)) can be ascertained in one of the 
ways described in ^2.41 

The scheme described in Sj2]can be used in a similar manner 
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Table 1. Intrinsic collision probabilities in 10 — 18 km -2 yr" 1 between the 
Astrid and objects from Table II of Dell'Oro & Paolicchi (1998) for com- 
parison of our results with those of that paper and with those of Bottke & 
Greenberg (1993). 



Object 


BG93 


DP98 


This paper 


1948 EA 


3.16 


3.19 


3.17 


Apollo 


3.58 


3.58 


3.58 


Adonis 


4.51 


4.52 


4.35 


1950 DA 


3.69 


3.76 


3.76 


Encke 


3.36 


3.42 


3.43 


Brorsen 


0.93 


0.94 


0.94 


Grigg-Mellish 


0.022 


0.022 


0.022 


Temple-Tuttle 


0.62 


0.62 


0.60 


Neujmin 


0.94 


0.93 


0.93 


Schaumasse 


1.13 


1.15 


1.15 


Pons Brooks 


0.041 


0.041 


0.041 



to work out the rate of catastrophic impacts onto planetesimals of 
size 7J max in cell (qi, Qi) from planetesimals in cell (g2, Q2) 

^cc(Anax, qi, Ql',12, Q2) = Af t ot(?2, Q2)Rcc(D max 5 Ql , Ql ', Q2 , Q2) 

For the 2-dimensional approximation the normalised collision rate 
is 

-Rcc(-Dmax, gi, Qi; 92, Q2) = A r y^g-(r)(gi,Qi) x 

r 

a v (r){q 2 ,Q 2 ){v 8 i ^ 12] )(r)i2S) 

where (v^^ 2 ]){ r ) ls me mean of the relative velocities to the 8/3 
power between planetesimals in the different cells at that radius 
(which can, e.g., be worked out using the Monte-Carlo method of 

The computation of -Rcc(fmax, Qi, Qi', 92, Q2) can be 
made more efficient noting that -R C c(Anax, gi, Qi; q2, Q2) = 
^cc(Anax, g2, Q2; 9i, Qi) and that some combinations of 
gi, Qi, g2, Q2 are either unphysical or non-overlapping, as well 
as by only including in the Monte Carlo calculation planetesimals 
where the orbits from the two cells overlap. A similar simplifica- 
tion to that used to derive equations J16M17t can also be employed 
so that collision rate is given by 

Rcc(D max ,qi,Qi;q 2 ,Q 2 ) = (K/2.1 x Kr 14 )M, 4/3 g- 13/3 x 

62, 92/(71). (29) 

Here-on we plot the function T(ei, 62,q2/qi) defined in this way 
rather than _R CC , noting that the situation described in H2.4\ cor- 
responds to T(e, e, 1) which must therefore be equal to 1.43 x 
HT^l - e) 13 / 3 5 2 d(e). 

3.2 Comparison with previous collision rate calculations 

The collision rates described in ^3.1l use the 2-dimensional approxi- 
mation, but can be readily modified for a 3-dimensional calculation 
in a manner similar to ^2.4.61 It is also possible to modify the as- 
sumptions to choose an inclination distribution that spans a narrow 
range 7 m j n to l ml «. In this case the collision rates we derive should 
be comparable with those of other authors who consider interac- 
tions between planetesimals on two orbits each of which is defined 
by a semimajor axis, eccentricity and inclination (with other angles 
assumed to be randomly distributed). 

To carry out such a comparison we also need to modify the 
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Figure 4. Contribution of different relative velocities to the intrinsic col- 
lision probability for two planetesimals, one with orbital elements a = 
3.42AU, e = 0.578, / = 0.435rad and the other with a = 1.59AU, 
e = 0.056, / = 0.466rad for comparison with Fig. 1 of Bottke et al. 
(1994) and Fig. 1 of Dell'Oro & Paolicchi (1998). 

•(27) 

calculation to derive the more typically quoted quantity of the in- 
trinsic collision probability, Pi, which is the probability of impact 
per unit time divided by r 2 where r = (D + D- lm )/2 (Wether- 
ill 1967). Our method included the r 2 factor from the outset (see 
equation O and integrated the collision rate over the size distri- 
bution capable of causing catastrophic impacts, whereas this was 
accounted for at a later stage by other methods. Using our notation, 
intrinsic probability in 10 -18 km -2 yr -1 is given by 

P = 0.030 ]T a(r, cf>)( qi , Qi)a v (r, 0)(g 2 , Qa)(tto[i,a]>(r, 0), (30) 

where the constant includes both a conversion between units and an 
extra factor it (since our calculation considers cross-sectional area 
rather than r 2 ). In Table[T]we compare the intrinsic collision proba- 
bilities we find using the Monte-Carlo approach with 200,000 plan- 
etesimals in each population with those of Dell'Oro & Paolicchi 
(1998) and Bottke & Greenberg (1993) showing excellent agree- 
ment between the methods. 

As noted by Bottke et al. (1994) (their point 1.2), to utilise 
intrinsic collision probabilities requires knowledge of the velocity 
probability distribution and it is not sufficient to assume a single 
mean relative velocity of collisions. This is automatically included 
in our method, since we keep track of the velocities encountered by 
a planetesimal at each location and these are used to work out the 
amount of cross-sectional area of impactors at that location that are 
able to cause a catastrophic collision. To illustrate that our Monte- 
Catio method computes an accurate velocity probability distribu- 
tion, here we reproduce Figure 1 of Bottke et al. (1994) (see also 
Fig. 1 of Dell'Oro & Paolicchi 1998) in our Figure|4] 

Similarly we find that we can reproduce Figs 5 and 6 of Krivov 
et al. (2005) (not shown here) noting that their A and t;i mp are given 
by 

A(a 1 ,e 1 ,a 2 ,e 2 ) = ^ a(r, cf>)(qi, Qi)a v (r, 0)(g 2 , Q2),(31) 

«imp = ^2a(r,4>)(qi,Qi)a v (r,<f))(q2,Q2) X 

{vret[i,2])(y, 0)/A(ai, ei, a 2 , e 2 ), (32) 
in our notation, although we should point out again that such a 
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mean collision velocity must be used with care as a range of ve- 
locities contribute to the collision rate (see Fig. [4}. 

In summary we conclude that, for the same assumptions, our 
method for calculating collision rates produces results that are in 
agreement with the most accurate methods available in the litera- 
ture, and that although this method was derived with the intention 
of studying high eccentricity orbits, it is also applicable regardless 
of the magnitude of the eccentricity. 

3.3 Implementation of evolution 

The factor fl cc (D mal , q%, Qi; 52, Q-i) (equation l28l> depends only 
on the way in which the parameter space is divided up. Thus an 
approach was implemented in which this factor was worked out 
ahead of time, and then used to work out the catastrophic collision 
rate for the largest planetesimals in cell (q, Q) from planetesimals 
in all other cells: 



M tot (q,Q y t+dt) = M t ot(g,0,t)[l-rft*i?cc(Ana X ,g,Q)].(34) 

The timestep dt is chosen so that some fraction, say 0.1%, of 
the mass is removed from the most rapidly evolving cell in that 
timestep, i.e., 

dt = 10~ 3 /max[R cc (D max ,q,Q)]. ( 35 ) 

3.4 Application to scattered disk-like distribution 

One of the simplest examples of this evolution is its application 
to a scattered disk-like distribution wherein planetesimals have a 
common pericentre distance, q, but a range of apocentre distances. 
Such a distribution might arise from planetesimals scattered by a 
planet on a circular orbit at a radius just inside q. 

Figure [5ji shows the function T(ei,e2, 1) for a population 
with pericentres in the range ( 1 ± 0. l)qi and apocentres in 30 loga- 
rithmically spaced bins in the range (1 — 50)gi. This illustrates how 
for most planetesimals in the population, their collision rate per 
unit mass of the colliding population is highest for collisions with 
the lowest eccentricity population. This is because the mass of a 
lower eccentricity colliding population is concentrated in a smaller 
volume resulting in higher densities and so higher collision rates. 
However, for the low eccentricity population it is noticeable that 
T(0, e2, 1) peaks at eccentricities of around 0.5, and that the col- 
lision rate per unit mass of the colliding population is lower with 
planetesimals on circular orbits. This is because the higher den- 
sity of the low eccentricity colliders is offset by the lower relative 
velocities. The reduction in relative velocities as eccentricity is de- 
creased was already noted in Figure |2j5, and its impact on collision 
rates evident in Figure |3p. 

To illustrate how such a population evolves through mutual 
collisions, we assume that the mass in orbits in the range Q to Q + 
dQ starts off oc Q~ a dQ with a total mass of Mdisko distributed 
among apocentres in the range Q/q = 1 — 50. This assumption 
results in an initial distribution of mass surface density that falls 
off approximately r Q , because most of the mass is 

concentrated at apocentre. Figure shows the evolution of mass 
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Figure 5. Collisional evolution of a scattered disk-like distribution, (a) The 
function T(ei, e2, 1) that through equation 1291 defines the rate of colli- 
sions for a planetesimal in population 1 interacting with planetesimals in 
population 2. (b) Evolution of the mass surface density distribution assum- 
ing an initial distribution of apocentres n(Q) oc Q~ Q ,wherea = l,2,3is 
shown with solid, dotted and dashed lines respectively. The resulting evolu- 
tion scales with q\ and Af^iskO as P er me axis labels, and 4 lines are shown 
for each initial apocentre distribution, at times of t' = 0, 10 4 , 10 6 , 10 8 , 
where real time in years is related to t' through equation (36). (c) Evolution 
of total disk mass for the distributions plotted in (b), including also a = — 1 
and 5. 



-Rcc(-Dmax, q, Q) = 22 ^2 Mtot (q2, Qi ) Rcc (Anax , q, Q; q 2 , Q 2 ).(33) 

12 Q2 

Using the assumptions about the evolution of the size distribu- 
tion described in £12.21 this collision rate can then be used to work 
out the mass which would be removed in a timestep dt, and the 
mass distribution stepped forward in time. 
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surface density for a = 1,2,3, where equation 1291 was used to 
normalise both radius and mass surface density, meaning that the 
times plotted correspond to real times of 

t = (2.1 x 10" 14 /-R')M- 4/3 (/ max /0.05)g 1 13/3 M d l s 1 k t' (36) 

years, where Mdisk is in M© and q\ in AU. 

The a — 1 evolution shows how the distribution tends to a 
flat distribution in mass surface density. This is because the low ec- 
centricity population is rapidly depleted, with mass becoming ever 
more concentrated in larger apocentre population that evolves rel- 
atively slowly; note from Fig. |2ji that a high eccentricity popula- 
tion would have a fairly flat mass surface density distribution. The 
a = 3 evolution differs from that of a = 1 in that a bimodal pop- 
ulation is formed with a low eccentricity population causing the 
peak in surface density at low r and a high eccentricity population 
causing the fiat surface density distribution. Because the mass starts 
off concentrated in the low eccentricity population, its evolution is 
unaffected by the material with large apocentres, and so it evolves 
due to mutual collisions. However, populations with higher eccen- 
tricities are rapidly depleted in collisions with the low eccentricity 
population, with the highest eccentricities persisting the longest. 
Eventually the low eccentricity population is depleted in mutual 
collisions so that the mass is concentrated in the highest eccentric- 
ities even for a ^ 3. 

The evolution of total disk mass is shown in Fig. |5p. From this 
it can be seen that the evolution is slower for shallower apocentre 
distributions. All distributions also tend to a mass evolution that 
falls off oc 1 /t, meaning that the amount of mass remaining at late 
times is independent of both initial mass and a. The amount of re- 
maining mass lies between two values: If the mass is in the lowest 
eccentricities (steep distributions at young ages) then that mass is 
given by that expected from a low eccentricity population at r = q. 
However, if the mass is in the highest eccentricities (all distribu- 
tions at late ages) then that mass is given by the mass expected for 
the highest eccentricities in the distribution; since e max = 0.96 in 
this simulation because Q/q ^ 50, this figure is in agreement with 
the expectations of Fig.[3p and[5^. Typically the evolution switches 
from having the mass in the lowest eccentricity population to the 
highest eccentricity population. 

Application of this can be readily seen for the case of the scat- 
tered disk in the Kuiper belt. For a pericentre at 30AU and an ini- 
tial mass of 0.1M® we see that the evolutionary timescale is very 
long, since real time is 25 x 10 6 t' (for K = 2.1 x 10 -14 and 
/max = 0.05), so that 4.5Gyr corresponds to t' = 200, mean- 
ing that collisional mass loss would be expected to be very small 
over the age of the Solar System, regardless of a. Had the scattered 
disk had 100 times more mass in the epoch prior to the Late Heavy 
Bombardment, real time would be 0.25 x 10*V, so that 800Myr 
corresponds to t' — 3200, and again we would not expect colli- 
sional evolution to be strong. However, it must also be noted that 
we have assumed here that planetesimals remain in their (q, Q) cell 
and are only removed by collisions. In the scattered disk dynamical 
processes, such as scattering by Neptune, can occur on timescales 
that are shorter than collisional mass loss, and so must be taken 
into account (e.g., Volk & Malhotra 2008), although the 2003 EL 6 i 
collisional family in the Kuiper belt may provide evidence of the 
role of collisions in the evolution of the scattered disk (Levison et 
al. 2008). 

Another application would be to a putative scattered disk 
around HD69830. For a pericentre at 1AU, just outside the out- 
ermost (known) planet (Lovis et al. 2006), we see that for an initial 
mass of 1M® real time is the same as t' (for K = 2.1 X 10 -14 and 



/max = 0.05), and so after ~ 2Gyr of evolution the maximum mass 
that can remain is independent of both initial mass and a, although 
in accordance with the earlier discussion it does depend on the max- 
imum eccentricity in the distribution. Thus we find that 10~ 3 M© 
can remain at 2Gyr if apocentres extend out to 50AU, with more 
mass remaining should higher apocentres be present (or lower val- 
ues of K be applicable). It is worth noting from Fig.[5p, however, 
that although we have increased the mass that passes through 1 AU, 
this does not necessarily increase the mass surface density at 1AU. 



4 EMISSION PROPERTIES OF ECCENTRIC 
PLANETESIMAL SWARMS 

It was shown in Sj2]and Sj3]how increasing planetesimal eccentricity 
leads to longer collisional lifetimes, and higher disk masses at late 
times in spite of collisional processing. Here we consider how the 
emission spectrum of a planetesimal swarm changes as its eccen- 
tricity is increased and how such swarms might appear in observa- 
tions of dust around nearby stars. 

4.1 Emission spectrum 

The emission spectrum from planetesimal belts comprised of dust 
with absorption (and emission) efficiencies Q a bs(A, D) can be cal- 
culated using: 

F v = 2.35xlO~ 11 d- 2 a tot J2a(r) Q abs (A,/J)x 

B u {\,T(D,r))a{D)dD, (37) 
T(D,r) = [(Qab s (D,A)) T ,/(Qab s ( J D,A)) T(I5 , r) ] 1/4 r bb (r),(38) 
T bb (r) = 278.3L°' 25 r -0 ' 5 , (39) 

where F u is in Jy, d is distance in pc, a to t is in AU 2 and scales 
with total mass according to equation {8]l for the assumptions of 
this paper about the size distribution, a(D)dD is the fraction of the 
total cross-sectional area in the size range D to D + dD and is oc 
D~ 15 dD for the distribution assumed here, Dmin is the minimum 
size in the distribution that is commonly assumed to be the size at 
which dust is blown out of the system by radiation pressure, and 
/)max is the size of the largest object, although for computational 
purposes this can be set at ~ lm as larger objects contribute little 
to the emission spectrum (e.g., Fig. 5 of Wyatt & Dent 2002), B v 
is in Jy sr _1 , (Q a bs)T means Q a bs averaged over a black body 
spectrum of temperature T, and L* is in Lq . 

Later in the paper we will use emission efficiencies calcu- 
lated using Mie theory (or another suitable approximation) along 
with the optical constants of different materials combined using a 
mixing theory in the manner described elsewhere (Li & Greenberg 
1998; Wyatt & Dent 2002). The grain model we will use assumes a 
core-mantle composed of silicates and organic refractory material, 
in the ratio qsi, which is incorporated into the grain with a porosity 
p; some fraction 9h 2 o of the vacuum is filled with water ice. How- 
ever, for heuristic purposes this section starts by assuming that dust 
created in the planetesimal belt both absorbs and emits light like a 
black body so that Q a bs = L 

For the disk comprised of planetesimals with common peri- 
centre and apocentre distances discussed in Sj2] Figure [6^ shows 
the emission spectrum for a range of eccentricities for a star with 
Li, = 11/0, d = lOpc, and for a planetesimal belt with q = 1AU 
and (Ttot = 10 _3 AU 2 . Although the absolute level of emission and 



Collisional evolution of eccentric planetesimal swarms 1 1 



range of wavelengths are dependent on these stellar and planetesi- 
mal belt properties, the shape of the emission spectrum would not 
be (e.g., wavelengths would scale oc LZ°' 2j q ' 5 ). This shows that 
increasing eccentricity leads both to emission over a wider range 
of wavelengths, due to the larger range of radii and so temperatures 
in the disk, and to a decrease in the level of emission for a given 
amount of cross-sectional area. 

To work out the emission spectra for these disks at late times, 
equation 1241 was used to get the mass, and equation [8] to get the 
cross-sectional area (see Fig.[6p). The minimum size in the distri- 
bution was assumed to be the size at which dust is blown-out by 
radiation pressure. Because the majority of collisions occur at peri- 
centre (see Figure[2j;), the blow-out size is that for which radiation 
pressure f3 — 0.5(1 — e), so that for black body grains 

D hl = 0.8(2700/p)(Z*/Af*)/(l - e), (40) 

in fim. The increase in blow-out size for high eccentricities is sig- 
nificant since it truncates the collisional cascade thus removing 
much of its cross-sectional area. Figure[6jj shows that the effect of 
increasing eccentricity of a planetesimal belt above 0.1 while keep- 
ing its pericentre distance constant is first to reduce the amount 
of emission present at late times at all wavelengths (compare the 
e = 0.1 and e = 0.6 lines). However, for high enough eccentric- 
ities the emission is increased above the e = 0.1 values, with the 
transition occurring at eccentricities of ~ 0.99 for short wavelength 
emission and ~ 0.7 for long wavelength emission. 

Although the quantitative conclusions about the absolute level 
of emission would be different if realistic particles had been as- 
sumed, these qualitative conclusions about how the emission spec- 
trum changes as eccentricity is changed would not, since a change 
in particle properties would affect the emission from all disks in a 
similar manner. To illustrate this, Fig.[6j: shows the emission spec- 
tra assuming that the total cross-sectional area is all in 1 /im par- 
ticles of silicate-organic refractory material (noting that this is not 
meant to be a physical model). The emission efficiencies of such 
particles drops rapidly at longer wavelengths (apart from close to 
the 10 and 18 pm silicate features) dramatically reducing the far-IR 
emission, even though the qualitative comparison of the behaviour 
as eccentricity is changed is unaffected (i.e., the eccentricity re- 
quired for an increase in short or long-wavelength emission is the 
same). 

4.2 Consequences of increased remaining mass 

While increasing the eccentricity does not have the effect of in- 
creasing the warm disk emission, it does increase the remaining 
mass, and that has two important implications. 

4.2.1 Blow-out population 

The first is that the collisional cascade is losing mass through 
radiation pressure blow-out at a rate M to t/icc(£ ) max ), which is 
A/i a te/£agc once the largest objects reach collisional equilibrium. 
This means that a higher remaining mass means a higher mass loss 
rate. Typically the short lifetime of blow-out grains means that they 
contribute little to the total cross-sectional area present in the disk. 
However, for an eccentric ring in which the majority of the mass 
loss occurs at pericentre, the surface density of the blow-out grains 
can exceed that of the collisional cascade. 

The following arguments give an estimate of the surface 
brightness of these grains. To simplify this calculation it is as- 
sumed that all collisions occur at pericentre where dust of size 
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Figure 6. Emission spectra of disks comprised of material all with the 
same pericentres, but with eccentricities of 0.1, 0.6, 0.9, 0.96, 0.99 shown 
with different linestyles. (a) Spectra for 10 _3 AU 2 of material with black- 
body emission properties and with pericentres at 1AU around a 1-Lq star 
at lOpc. (b) As for (a), but the total cross-sectional area has been scaled 
to the maximum mass that can remain after lGyr of processing (assuming 
K = 2.1 X 10 -14 and / max = 0.05), also taking into account the change 
in blow-out radius as eccentricity increases, (c) As for (b), but assuming 
the cross-sectional area is in 1 fim grains of a mixture of silicate-organic 
refractory material. 
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Di OBB (in /im) is released at the orbital velocity of v\ OBB = 
27r x /M*(l + e)/g AU yr -1 . Further assuming that this dust has 
(3=1 means that the dust undergoes no acceleration so that 

n oss = uioss\A - (q/r) 2 . (41) 

The mass loss rate, when converted to a rate of loss of cross- 
sectional area as dust of size Dioss and divided by the cross- 
sectional area present in the collisional cascade gives 

o-ioss/o-totcc = 3.16 x 10 4 D°f D^D^t^, (42) 
which results in a distribution of cross-sectional area for r > q of 

doi oss /o totcc /dr = 4500L°' 5 A/~ 1 Z)^xAos^ag 1 c'3 1 ' 5 x 

(l + eJ-^fl-er^Il-Cg/r) 2 ]!^) 

This can be compared with the distribution of cross-sectional 
area in the collisional cascade (Fig. |2^), which the analytical re- 
sults (eq. [T8l show has a minimum at f = 1 — e 2 of a/df = 
n~ 1 \/e- 2 — 1. Thus the blow-out population becomes more im- 
portant relative to the collisional cascade as eccentricity increases, 
and is also more dominant in populations with larger pericentre dis- 
tances. Specifically, the cross-sectional area in the blow-out popu- 
lation can be higher than that of the collisional cascade at some 
radius when 

e > 1 - 160L? 25 M-° V- 75 ^SA; a ?C 5 (44) 

for these assumptions. For the A/* = = q = 1 system at lGyr 
considered earlier, further assuming that D max = 2000km and 
Di oss = 1/im, this means that the blow-out population becomes 
important when e > 0.97, noting that D\ OSB <C Dbi = 30/im in 
this example. 

4.2.2 Cometary sublimation 

The conclusions of ^4.2. II are independent of the mechanism pro- 
ducing the blow-out grains, which need not be through collisions. 
Rather the large remaining mass might provide a reservoir that, 
through collisions, replenishes smaller planetesimals that release 
dust through comet-like sublimation. In such a case the above cal- 
culation may overestimate the amount of cross-sectional area in the 
collisional cascade meaning that the cascade could be extremely 
faint, because sublimation timescales (which are oc 73 if mass loss 
rate is proportional to surface area) can be shorter than collision 
timescales (which from eq. [4] are oc D 0,5 ) for the small particles 
in the cascade, thus reducing their number so that D m - m should 
be significantly higher than D b i in equation l[8j. However, this 
would only be the case if sublimation processes just produce non- 
sublimating (e.g., silicate) grains that are below the blow-out limit. 
Large non-sublimating grains produced above the blow-out limit 
would increase the cross-sectional area in the collisional cascade. 
Since Solar System comets are seen to release significant quantities 
of non-icy mm-cm-sized grains (e.g., Reach et al. 2009), we con- 
sider that further work is required before the impact of sublimation 
on the steady-state size distribution of a collisional cascade is fully 
understood, but note that sublimation remains a viable mechanism 
for feeding the blow-out population discussed in £14.2.11 

4.2.3 Frequency of massive collisions 

To estimate the frequency of massive collisions for a system evolv- 
ing in collisional equilibrium, for which the collisional lifetime of 
the largest objects is the age of the system, we rewrite equation 24 



of Wyatt et al. (2007a) to find that the collision rate for objects 
larger than size D pb is given by 

R CC {D > D pb ) = (Mi a te/6M pb )t-e, (45) 

where a distribution with q — 11/6 was assumed. This means that, 
for a given parent body size, a higher eccentricity ring results in a 
higher remaining mass and so more frequent collisions. 

However, a higher eccentricity also means that larger parent 
bodies are required to reproduce the same fractional luminosity. To 
estimate the minimum mass of a parent body, D p h, that would be 
required to be destroyed to reproduce an observed fractional lumi- 
nosity of fobs, it is assumed (optimistically) that a collision turns 
all of an object's mass into particles of size Dm which are then re- 
distributed around a ring with the same eccentricity as the parent 
object (i.e., ignoring the increase in eccentricity due to radiation 
pressure). The fractional luminosity from an eccentric ring can be 
calculated by integrating equation dl8t to be 

jobs = cT tot /(4nq 2 {l - erVl - e 2 )- (46) 

Thus to reproduce / obs requires a parent body mass 

M pb = 3.14 x 10" 8 ^bi/ob s g 2 (l - e)" 2 (l - e 2 ) ' 5 (47) 

in Mq, noting that D^i is also a function of eccentricity. This read- 
ily gives the maximum rate of collisions that could reproduce a 
given fractional luminosity in a single event; e.g., if a mass of lAf® 
remains at lGyr in a belt with a pericentre of 1AU and eccentricity 
of 0.99, then events that might potentially produce / bs = 10 -3 in 
1/im dust occur 1.4 times per Myr. 

On the other hand a higher eccentricity means that the dust 
is not depleted in collisions so rapidly. The collision rate of same 
sized particles can be worked out using equation © with f cc (D) = 
4 (Wyatt et al. 1999). Using the analytical approximation of £ |2.4.2| 
we find that 

Rcc(D bl ) = 1.3a tot M?- V 3 ' 5 (l - e) 3 ' 5 (l - e 2 )- a5 / n "L, (48) 

which is true regardless of particle size (as long as such particles 
contain the majority of o- to t), and gives a collision rate to reproduce 
UsOfWUsM^I-^a- 1 - 5 . 

Combining these results in the same manner as Wyatt et al. 
(2007a), to estimate the fraction of time that collisions are expected 
to result in dust above a given level of Jobs, we find that 

P(f > Us) = 0.33 x 10 6 M lato t- 1 o /- b 2 AC°' 5 / max x 

p-^biV^a+er - 5 . (49) 

For example, if a mass of 1M® remains at lGyr in a belt with a 
pericentre of 1AU and eccentricity of 0.99, then events that might 
potentially produce / Q b s = 10~ 3 in 1/im dust of density 2700 kg 
m -3 with inclinations up to 0.05, would occur 0.4% of the time. 
This fraction is insensitive to eccentricity except that a high value 
is required to permit 1M® of material to remain so late despite col- 
lisional processing, and would also affect the validity of assuming 
the dust is 1 in size. 

4.3 Application to r/ Corvi 

The 1.3 Gyr F2V star rj Corvi at 18.2pc exhibits excess emission 
from circumstellar dust at wavelengths from a few /im up to sub- 
mm. The sub-mm emission has been resolved in imaging at both 
450 and 850/im with SCUBA, and modelling shows this dataset 
can be explained by emission from a ring of 150 ± 20 AU radius 
inclined at 45 ± 25° to our line-of-sight (Wyatt et al. 2005). The 
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Figure 7. Emission spectrum of r\ Corvi modelled using a single planetes- 
imal population with pericentres at 0.75AU and apocentres at 150AU (see 
text for details). The dashed, dotted and solid lines correspond to the con- 
tribution from the collisional cascade, from blow-out grains, and the total 
emission spectrum, respectively. The diagonal solid line is the stellar spec- 
trum, the grey dotted line is the IRS spectrum after the photospheric contri- 
bution has been subtracted, asterisks and squares are the IRAS and Spitzer 
photometric fluxes, respectively, after photosphere subtraction (and colour 
correction where necessary). 



region is seen to be centrally cleared, but the excess mid-IR emis- 
sion cannot originate in the 150AU ring and has been resolved to 
lie somewhere in the 0. 16-3.0AU region (Smith et al. 2008; Smith, 
Wyatt & Haniff 2009b) with a temperature inferred from the emis- 
sion spectrum that places it at ~ 1.7AU. It remains debatable 
whether there is any emission in the 3-100AU region (Chen et al. 
2006; Smith et al. 2008). Regardless the warm emission at < 3AU 
is expected to be transient if confined to a ring at that location (Wy- 
att et al. 2007a), as discussed in the introduction. 

Here we consider the possibility that both the imaging and 
spectroscopic constraints can be explained by a model with a single 
planetesimal population. Figure [7] shows the emission spectrum of 
a model with planetesimals of density 1480 kg m~ 3 with pericen- 
tres at 0.75AU, apocentres at 150AU (i.e., eccentricity 0.99), and 
inclinations up to 0.05, after 1.3Gyr of evolution. The planetesimal 
collisional parameters are assumed to be D max = 2000km with 
strength — 3 x 10 5 J kg -1 , leading to the conclusion that 
of material remains at the current epoch. Dust is created in the colli- 
sional cascade down to a size of Dbi = 480 before removal by 
radiation pressure, and the emission properties of the grains were 
calculated assuming a core-mantle composed of 30% amorphous 
silicate and 70% organic refractory material with porosity of 0.4 
and water ice filling 20% of the gaps; the emission spectrum of the 
collisional cascade is shown with a dashed line, and is largely in- 
sensitive to the assumed composition as bound grains act like black 
bodies. The population of grains that are being removed by radia- 
tion pressure are assumed to be all of 5/im in diameter, and their 
emission properties were calculated assuming the same grain prop- 
erties as the collisional cascade; their emission spectrum is shown 
with a dotted line. 

The model provides an excellent fit to the observed shape of 
the emission spectrum, noting that it is not our intention to pro- 
vide an exact fit to the emission features for which a more de- 
tailed compositional model would be required; i.e., the composi- 
tion should not be regarded to be constrained in this model, but in 
providing a fit to the spectrum it does allow the physical param- 



eters of the model to be self-consistently derived. In the model, 
the mid-IR emission is dominated by grains being removed by ra- 
diation pressure, which contribute 91% and 84% of the 11.6/im 
and 18.7/im emission, respectively. Such emission is centrally con- 
centrated (e.g., 68% and 51% of the emission in these wavebands 
comes from inside 3.5AU). To assess whether the model is con- 
sistent with available mid-IR imaging constraints we first show 
that the model fluxes in an aperture of radius 0.5"(285mJy and 
263mJy at these wavebands) are consistent with those observed 
(330 ± 184mJy and 309 ± 79mJy, Smith et al. 2008). In the model, 
the aperture used to estimate the background level, 0.5-1.0", does 
include disk emission (28mJy and 60mJy in these wavebands), but 
at a level consistent with its non-detection in Smith et al. (2008) 
for which > 150 and > 140mJy would have been necessary for 
a > 3cr detection. To consider whether there is an aperture which 
could confidently detect the mid-IR emission from such a model 
we plotted the model fluxes on Fig. 7 of Smith et al. (2008), noting 
that the central concentration of the mid-IR emission means that we 
only need to consider the emission close to the pericentre. Thus we 
approximated the mid-IR emission from the model as originating in 
two consecutive annuli each with a uniform surface brightness and 
spanning radii 3 ± 3AU and 9 ± 3AU; from the model we find that 
the inner annulus contains 262mJy and 239mJy, and the outer an- 
nulus contains 27mJy and 53mJy at these wavebands. Both annuli 
are below the threshold at which extended emission should have 
been detected, consistent with the observations, but the outer annu- 
lus is at a level where extension would be detected in observations 
twice as deep at both wavebands. The lack of material < 0.75AU 
is consistent with visibilities measured by MIDI that suggest the 
emission is completely resolved, although the model would predict 
small changes in visibility with baseline for a 40mas ring that might 
be detectable in more sensitive observations (see Fig. 9 of Smith et 
al. 2009b). The low emission efficiencies of the blow-out popula- 
tion means that it contributes little to the longer wavelength fluxes, 
and the sub-mm emission is dominated by the collisional cascade. 
The majority of the sub-mm emission comes from material near 
apocentre (92% of the 450 /im emission comes from 100 — 150 AU) 
and so we expect the model to fit the sub-mm imaging constraints 
(Wyatt et al. 2005). 

Although we have devised a model that explains the obser- 
vations, this does not necessarily make it plausible. Nevertheless, 
most of the model parameters are reasonable. The planetesimal 
strength required is above that assumed for the population models 
of Wyatt et al. (2007b), but is within the range of that expected for 
2000km planetesimals (Benz & Asphaug 1999). The assumption 
that blow-out particles are 100 times smaller than the blow-out limit 
could indicate that larger unbound particles readily disintegrate into 
smaller fragments (noting that they are produced at collisional ve- 
locities up to 80 km s^ 1 ), but may perhaps be circumvented with 
a different choice of composition. It is also worth noting that this 
model predicts that we have a 1 :2300 chance of witnessing a col- 
lision capable of reproducing the observed fractional luminosity of 
0.5 x 10~ 3 ; i.e., it is possible that the mid-infrared emission is en- 
hanced due to stochastic collisions. 

The origin of the planetesimals is, however, a concern. High 
eccentricity orbits normally imply scattering by a planet, with the 
apocentre or pericentre coinciding with the orbit of the planet. It 
is unclear why both pericentre and apocentre distances would be 
fixed here. One possibility is that a primordial population that in- 
cluded a wider range of eccentricities and semimajor axes has been 
depleted by collisional processing leaving just the high eccentric- 
ity remnant (see $3.4\ . The primordial population could then be an 



14 M. C. Wyattetal. 



extended scattered disk created by inward planet migration (Payne 
et al. 2009), with a planet lying just interior to the pericentre, and 
the apocentre corresponding to the highest eccentricity attained in 
the population. In this case the current mass is a small fraction of 
the primordial mass; e.g., for n(Q) oc Q~ a where a = 1 or 2 
we find the primordial mass to be 7 or 200 times (respectively) the 
5M® inferred at present (assuming the observed apocentres are in 
the range 90-210AU), and so an initially flat apocentre distribution 
is required for a realistic starting population. Alternatively this pop- 
ulation could have been scattered in by a more distant planet which 
was orbiting close to the apocentre, although in this scenario the 
lack of material with pericentres in the 3-100AU range is a prob- 
lem, as collisional processing preferentially removes material with 
the lowest pericentres in a distribution with a common apocentre. 
Certainly the dynamics of creating populations of extremely high 
eccentricity must be explored before we can be confident that this 
is a viable model. It is also important that the dynamics that creates 
the high planetesimal eccentricities, such as scattering by a planet, 
does not significantly deplete the planetesimal population on longer 
timescales. This could be achieved if the planetesimals become de- 
tached from the planet, perhaps because that planet migrates or is 
scattered to put it out of reach of the planetesimals. 

Regardless of the origin of the planetesimals, it seems that the 
0.75-150AU region would have to be clear of planets for this model 
to be viable. Current limits from radial velocity studies only rule out 
companions down to 2.1 Jupiter mass out to 0.48AU (Lagrange et 
al. 2009). As noted above, another testable prediction of the model 
is that there should be thermal emission from the 3-100AU region 
which may be detectable as extended emission. 



4.4 Application to HD69830 

The 2Gyr K0V star HD69830 at 12.6pc exhibits excess emission 
from circumstellar dust at mid-infrared wavelengths, including sig- 
nificant emission features (Beichman et al. 2005). No excess emis- 
sion is seen at 70/im (Trilling et al. 2008) or at 850/im (Sheret, Dent 
& Wyatt 2004; Matthews, Kalas & Wyatt 2007). Three Neptune- 
mass planets are known within 0.7AU (Lovis et al. 2006) that are 
thought to have migrated in from 3-8AU (Alibert et al. 2006). 
Mid-infrared studies constrain the dust to 0.05-2.4AU (Smith et al. 
2009b), consistent with the temperature inferred from the emission 
spectrum that places the dust at 1AU (Lisse et al. 2007); at such a 
location the dust is inferred to be transient (Wyatt et al. 2007a). 

The lack of cold emission rules out several models for the ori- 
gin of the dust, including the eccentric planetesimal belt model as 
applied to r/ Corvi above. Although an eccentric ring alleviates con- 
cerns about the longevity of the disk feeding the dust at 1AU, the 
hot dust is always accompanied by cold emission from collisional 
cascade material at apocentre (Fig. |6p). One way around this is to 
postulate that there is little cross-sectional area in the collisional 
cascade, e.g. because small grains have been removed by sublima- 
tion ( i)4.2.2| >. Figure [8] shows a model which fits the emission spec- 
trum by only including the blow-out component comprised of solid 
grains 1/im in diameter composed of 1/3 amorphous silicate and 
2/3 organic refractory material; note again that the model is used 
only to fit the overall shape of the spectrum, and does not claim 
the level of detail required to constrain the composition. Dust in 
this model extends out from the 1.2AU pericentre of a ring with 



2 Collision rates for material with common apocentres, Qi, would be given 
by equation (29) replacing q^ 13/3 with Q^ 13/3 [(1 + ei)/(l - ei)] 13 / 3 . 
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Figure 8. Emission spectrum of HD69830 modelled using a single planetes- 
imal population with pericentres at 1 .2AU and eccentricity 0.99 (see text for 
details). The solid line shows the contribution from blow-out grains; for this 
to be a viable model the emission from the collisional cascade component 
has been suppressed, e.g., due to the sublimation of its smallest grains. The 
diagonal solid line is the stellar spectrum, the grey dotted line is the IRS 
spectrum after the photospheric contribution has been subtracted, asterisks 
and squares are the IRAS and Spitzer photometric fluxes, respectively, after 
photosphere subtraction (and colour correction where necessary). 



eccentricity 0.99; the spectrum is insensitive to eccentricity, how- 
ever this does have a moderate impact on the inferred mass loss 
rate through equation d4 1 1 . As well as nearly fitting the 70/im 
upper limit, this model also meets the mid-infrared imaging con- 
straints, since it predicts the 18.7/im flux within a 0.5"aperture to 
be 412mJy (compared with 377 ± 46mJy found by Smith et al. 
2009b), with just 30mJy coming from the 2-4" region. As for ^431 
the model was analysed assuming the emission to originate in two 
annuli at 2 ± 2AU and 6 ± 2AU with disk fluxes of 121 and 29mJy, 
respectively. Comparison with Fig. 7 of Smith et al. (2009b) shows 
that both annuli are below the threshold at which extended emission 
should have been detected in these observations, but that extension 
would have been detected from the outer annulus in observations 
twice as deep, a conclusion which holds regardless of the origin of 
the dust as long as it is assumed to be in the process of radiation 
pressure blow-out. Without material < 1.2AU the model emission 
would be expected to be completely resolved on MIDI baselines in 
line with observations (Smith et al. 2009b). 

This model gives a mass loss rate of 0.08M® Myr -1 , which 
is higher than the 5 — 60 x 10 _6 M^/Myr quoted in Beichman 
et al. (2005) and Wyatt et al. (2007a), since those papers assumed 
mass loss due to collisions rather than radiation pressure. Assum- 
ing this mass loss comes from a parent population in collisional 
equilibrium, this implies a parent population of I6OM0, similar to 
the 90 — 900M® calculated in Beichman et al. (2005) based on an 
extrapolation from Hale-Bopp mass loss rates. It seems that such a 
high remaining mass is prohibitive after 2Gyr of processing, since 
this can only be achieved for e > 0.998 for the fiducial value of 
K = 2 x 10" 14 or e > 0.993 for the slightly higher K-value used 
for r\ Corvi modelling (i ]4.3t . Even if it was possible to implant 
such a large mass at high eccentricities, there would still be the 
concern raised in ^4.2. 21 about whether comet sublimation would 
also produce bound non- sublimating grains that would increase the 
cold emission from the collisional cascade. Thus a steady-state ex- 
planation is not favoured. 

An alternative model invokes a recent collision as the origin 
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of the dust. In an eccentric ring such events can occur frequently 
enough for this model to be viable ( i]4.2.3t . A planetesimal belt 
mass of lM ffi results in P(f > 0.2 x 1(T 3 ) = 0.06 assuming 
that collisions result in 1/im grains. Since the observed grains also 
seem to be small enough to be blown out by radiation pressure, the 
probability of witnessing such an event is in fact lower than that 
derived in ^4.2.31 and quoted above: a short grain lifetime either 
implies the collision occurred very recently (i.e., within the blow- 
out timescale of several years), or that the dust is the product of 
ongoing secondary collisions or sublimation amongst the debris of 
an older more massive (and hence rarer) collision. The persistence 
of 24 /j,m excess over 24 year timescales between the epochs of 
IRAS and Spitzer (Lisse et al. 2007) rules out a very recent ori- 
gin. This leaves ongoing secondary collisions or sublimation in a 
debris belt which would have to contain > 2 x 1O _6 M0 to sus- 
tain a mass loss of 0.08M® Myr -1 over 24 years. Equation d45t 
shows that such debris belts might only be recreated in collisions 
every 24,000 years for a planetesimal belt containing 1 Af® at 2Gyr. 
In other words, we have a 1:1000 chance of witnessing the after- 
math of such an event, a conclusion which holds even if the mass 
loss had been assumed to persist over longer timescales (since the 
increased debris belt mass required would have been recreated in 
collisions correspondingly less frequently). Since 1M® is not an 
implausible remnant mass for an eccentric planetesimal belt (see 
$3.41 , we consider secondary collisions to be a possible explana- 
tion for HD69830, even if some process like sublimation must be 
invoked to deplete small grains and render the collisional cascade 
of the parent planetesimal belt non-detectable. 



5 CONCLUSION 

This paper considers collisional processes in a population of plan- 
etesimals with high eccentricities. Collision rates are derived both 
analytically and using Monte-Carlo simulations. It was found that 
eccentricity has a significant effect on collision rates, and that the 
amount of mass that can remain in a planetesimal belt at late times 
can be significantly increased. 

The emission properties of eccentric planetesimal belts were 
presented, and it was found that radiation pressure causes eccentric 
rings to be deficient in small particles. Thus, despite the increased 
mass of a high eccentricity planetesimal belt at late times, extreme 
eccentricities of > 0.99 are required to enhance the emission level 
above that expected from a low eccentricity belt. However, the high 
mass loss rate of extreme eccentricity planetesimal belts can cause 
the wind of blow-out particles that extends outward from the peri- 
centre to be detectable. The high frequency of massive collisions 
in such belts can also make it likely for us to be witnessing dust 
produced in such collisions. 

Application of this model to the rj Corvi debris disk showed 
that all available observations can be explained by an extreme ec- 
centricity (e = 0.99) planetesimal belt of mass 5M®, circumvent- 
ing the conclusion that the hot dust at 1.7AU must be transient. 
Despite this success, the dynamical challenges of creating such a 
massive extreme eccentricity population would need to be over- 
come before this model can be considered viable. Observational 
tests are suggested including the presence of resolvable emission 
(and absence of planets) in the 3-100AU region. Application to 
HD69830 is complicated by the lack of far-infrared emission. It 
may be possible to explain this non-detection by the removal of 
small dust from the collisional cascade by sublimation, in which 
case the mid-infrared emission may be plausibly explained by the 



ongoing destruction of a debris belt produced in a recent collision 
in an eccentric planetesimal belt. 

Although the majority of the discussion focuses on single 
eccentricity populations, the results can also be applied to pop- 
ulations with a range of semimajor axes and eccentricities. This 
was demonstrated by application to scattered disk-like populations 
where it was found that, in the absence of other dynamical pro- 
cesses, rapid collisional erosion of the low eccentricity populations 
would be expected to result in a single high eccentricity popula- 
tion. Since the known high eccentricity planetesimal populations 
are produced in interactions with planets, and so may be contin- 
uing to undergo dynamical evolution on timescales shorter than 
collisional timescales, it is noted that dynamical interactions may 
continue to play a defining role in the long term evolution of high 
eccentricity populations, and that the collisional evolution scheme 
presented here could be readily incorporated into N-body simula- 
tions of planet-planetesimal interactions to derive simultaneously 
the collisional and dynamical evolution of a scattered planetesi- 
mal population. A further extension of the model would include 
a prescription for planetesimal strength as a function of size which 
would lead to a departure from the single phase size distribution 
assumed here (e.g., Lohne et al. 2008). 

The results of this study would be applicable wherever non- 
negligible planetesimal eccentricities are found. Thus, other poten- 
tial applications include the Solar System's comet and NEO pop- 
ulations, particularly in the early phases when these populations 
would have been more massive and hence collisional processing 
more important (Booth et al. 2009), the outcome of planet forma- 
tion models (e.g., Payne et al. 2009), and systems where eccentric 
planetesimals may be implicated such as the origin of dust around 
White Dwarfs (Farihi, Jura & Zuckerman 2009) and of the hottest 
dust population of debris disks like Vega (Absil et al. 2006). 
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